Working with Spatial Data in terra and sf

sf is a new package and apepars to be the leading contender to replace rgdal and sp. sf reads and writes shapefiles and also has a large number of vector transforms available (see resources at the bottom of the page).

In the code below, you'll notice that R orders geographic coordinate values as "longitude, latitude" rather than lat,lon. This is common in spatial programming as we order everything as x,y and typically we treat longitude as an x and latitude as a y when working with geographic data. This also alighns with Easting, Northing.

The rspatial.org website is the best resource I have found for this material but it does not have the best tutorials. There are a lot of other websites out there but be careful as many of them are out of date.

Working with Vector Data

Loading a Shapefile

The code below will read a shapefile into a dataframe.

# Load libraries
library(terra)
library(sf)

# set the working directory so we don't have to specify a full path each 
# time we load/save a file
setwd("C:/Users/Jim/Desktop/GSP 570/Lab 4 GLM")

# load the shapefile into a SpaVector object
TerraVectorObject=vect("test.shp")

# use the standard plot function to show the spatial data
plot(TerraVectorObject)

# get the number of features in the shapefile
NumFeatures=nrow(TerraVectorObject)

# get one column of attribute values 
ColumnOfAttributeValues=TerraVectorObject$measyear

Because a SpatVector object is similar to a dataframe, the dataframe functions such as nrow() and ncol() work with the object.

Saving a Shapefile

writeVector(TheSpatVectorObject, "TestWriteTerra.shp")

Loading a CSV with Points

TheSpatVectorObject=vect(DataTableOfPoints,geom=c("lon", "lat"),crs=CRSDef)
plot(TheSpatVectorObject)

Saving a SpatVector object to a CSV

 

Accessing the Spatial Data

The SpatVector object handles points, polylines, and polygons. We'll mostly be using points. We can access the spatial data using the 'geom()' function. This will return a dataframe that contains the x (lon or easting) and y (lat or northing) values in columns 3 and 4. Then, we can access each feature's coordinate by indexing into the dataframe as below.

TheGeom=geom(TerraVectorObject) # provides a matrix of values

# get the coorindate for the second feature in the file
Lon=TheGeom[2,3] # column 3 are the lons
Lat=TheGeom[2,4] # column 4 contains the lats

 

Cropping a Raster

Cropping a raster is realively easy once you are used to working with extent rectangles. The "ext()" function will return the bounds of a raster in spatial coordinates (lon/lat or easting/northing) in an array of 4 values. The values are ordered as:

1. xmin
2. xmax
3. ymin
4. ymax

The code below gets a raster's extent and then modifies it to crop it in the x and y directions. You can plot a raster to see it's extent to generally crop to an area or use a GIS application to find the coordinates for a crop more precisely.

#get the extent object from an existing SpatRaster object
TheExtent=ext(TheRaster)
TheExtent # print the extent

#constrain it in the x direction by chainging the longitudes
#constrain it in x direction
TheExtent[1]=-124.75 # west 
TheExtent[2]=-114 # east
TheExtent[3]=30 # south
TheExtent[4]=45 # north
CroppedRaster=crop(TheRaster, TheExtent)
plot(CroppedRaster)

Converting a raster to a dataframe

Converting rasters to dataframes and back again is important in R as most modeling packages expect a dataframe rather than a raster. Passing "xy=TRUE" will include the xy coordinate values in the dataframe.

DataFameWithRasterData=as.data.frame(ProjectedRaster,xy=TRUE)
DataFameWithRasterData

Projecting Raster and Vector Objects

Projecting SpatVector and SpatRaster objects is easy once you have the proj string to define the projection.

ProjectedRaster=project(CroppedRaster, "+proj=utm +zone=10")
plot(ProjectedRaster)

 

Coordinate Reference Systems

You can use Proj strings to create Coordinate Reference Systems (crs) (a.k.a. Spatial Reference Systems). The example below creates an unproject (geographic) crs and then create a vector dataset that uses that crs.

CRSDef <- "+proj=longlat +datum=WGS84" # proj string

pts <- vect(lonlat, crs=crdref)

Other proj strings include:

+proj=utm +zone=10 # UTM Zone 10 North (North is the default)
+proj=utm +zone=59 +south # UTM Zone 59 South

# Robinson projection method. Try changing the longitude of origin
+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m # robinson in meters

The state plane CRSs have different projection methods and parameters. This information is available online at a number of websites (see below). As an example, State Plane California uses the Lambert Conic Conformal project method. Which has the following parameters (maptools.org).

+proj=lcc   
+lat_1=Latitude of first standard parallel             
+lat_2=Latitude of second standard parallel               
+lat_0=Latitude of false origin (i.e. latitude of origin)                 
+lon_0=Longitude of false origin (i.e. longitude of origin)
+x_0=False Origin Easting                
+y_0=False Origin Northing

The proj string for NAD 83 California State Plane Zone I meters would then be (EPSG.io):

+proj=lcc +lat_1=41.6666666666667 +lat_2=40 +lat_0=39.3333333333333 
+lon_0=-122 +x_0=6561666.667 +y_0=1640416.667

You may see the "+no_defs" parameter in some examples. This is an old parameter that prevented proj from reading a defualts file. This feature no longer exists in proj so "no_defs" it not needed.

You also may see "+towgs84=0,0,0" or something similar. This defines a shift between datums but if the values are 0 then it is not needed.

Additional Resources

Simple Features for R - website for SF with great quick reference.

Spatial Data Science - Vector Data - great detailed explaination of tarra.